The formulas and their derivation are from
>function f(x,y) ...
th1=y[1]; th2=y[2]; dth=th1-th2; c=cos(th1-th2);
l=1; g=10; m=1; ml2=m*l^2;
pth1=y[3]; pth2=y[4];
Dth1=6/ml2*(2*pth1-3*c*pth2)/(16-9*c^2);
Dth2=6/ml2*(8*pth2-3*c*pth1)/(16-9*c^2);
return [Dth1,Dth2,
-ml2/2*(Dth1*Dth2*sin(th1-th2)+3*g/l*sin(th1)),
-ml2/2*(-Dth1*Dth2*sin(th1-th2)+g/l*sin(th2))]
endfunction
>t=0:0.01:50;
>s=ode("f",t,[0,0,4,4]); s=s[1:2]-pi/2;
The next plot shows all positions reached from the end of the second axis.
>plot2d(cos(s[1])+cos(s[2]),sin(s[1])+sin(s[2])):
>t=0:0.05:50;
>s=ode("f",t,[0,0,4,4]); s=s[1:2]-pi/2;
We animate this.
>function anim (th1,th2) ...
setplot(-2,2,-2,2);
loop 1 to cols(th1)
clg;
a=th1[#]; b=th2[#];
lw=linewidth(3);
hold on;
plot([0,cos(a),(cos(a)+cos(b))],
[0,sin(a),(sin(a)+sin(b))]);
markerstyle("o#");
mark(0,0);
hold off;
linewidth(lw);
wait(0);
if testkey() then break; endif;
end;
endfunction
Press any key to sto the animation.
>anim(s[1],s[2]):
I have uploaded the animation to YouTube.